Entropy and density of states from isoenergetic nonequilibrium processes 
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Two identities in statistical mechanics involving entropy differences (or ratios of density of states) 
at constant energy are derived. The first provides a nontrivial extension of the Jarzynski equality 
to the microcanonical ensemble [C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997)], which can be 
seen as a "fast-switching" version of the adiabatic switching method for computing entropies [M. 
Watanabe, W. P. Reinhardt, Phys. Rev. Lett. 65, 3301 (1990)]. The second is a thermodynamic 
integration formula analogous to a well-known expression for free energies, and follows after taking 
the quasistatic limit of the first. Both identities can be conveniently used in conjunction with a 
scaling relation (herein derived) that allows one to extrapolate measurements taken at a single 
energy to a wide range of energy values. Practical aspects of these identities in the context of 
numerical simulations are discussed. 
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I. INTRODUCTION 

In the ubiquitous studies of phase stability, phase tran- 
sitions and reaction directions, most properties of in- 
terest are embodied in the free energies and entropies 
corresponding to the desired phases or reaction states. 
Not surprisingly, a great amount of intellectual effort has 
been directed towards finding ever more efficient compu- 
tational means of estimating these quantities . Within 
the wider domain of equilibrium thermostatistics, a sim- 
ilar role is played by the so-called density of states, the 
knowledge of which generates not only free energies and 
entropies, but also heat capacities and the very partition 
function. Here, too, one finds a great variety of computa- 
tional techniques specifically designed for its estimation 
(see, e.g., and references therein). 

Traditionally, the aforementioned computational 
methods are inherently nondynamical, i.e. the rele- 
vant data is obtained from either a single equilibrium 
macrostate, or a series of such states. An exception is 
provided by the adiabatic switching method of Watanabe 
and Reinhardt which allows one to recover entropy 
differences from a single dynamical trajectory connecting 
the states of interest. In principle, this method yields an 
exact estimate of the entropies provided the "switching 
process" is infinitely slow; in practice, one has to cope 
with the systematic errors that arise due to the unavoid- 
able finite switching time 0] (see also Q for further 
aspects of this method). In the context of free energy 
differences, this situation has substantially changed after 
an identity connecting free energy differences at constant 
temperature and nonequilibrium processes was derived 
by Jarzynski (^|6|], viz. 



-AF(T)/T 



(1) 
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where AF(T) = Fb{T) - Fa{T) is the (Helmholtz) free 
energy difference between the states of interest A and B, 
the angle brackets denote an average over all realizations 
of a predefined switching process connecting the states A 
and B, and W is the net amount of work performed on 
the system during each realization (throughout this pa- 
per, the units will be chosen so that the Boltzmann con- 
stant fcs equals unity) . Unlike the method of Ref. |1| , the 
above equality does not depend upon the rate at which 
the dynamical switching process is performed, and hence 
is free from the systematic errors associated with finite 
switching times. Another welcome computational aspect 
of this identity is that it is trivially parallelizable - each 
realization of the process can be performed separately. 
Because of these features, the Jarzynski equality (JE) is 
becoming increasingly popular in the simulation commu- 
nity, and various reviews are already available (see e.g. 

9]). 

Although its generalization to free energies other than 
the Helmholtz one is straightforward [3, [l3 > the analog 
of the JE for entropies cannot be obtained by such a 
direct route. In this case, one is intuitively tempted to 
adapt the derivation of Ref. 5] by simply replacing the 
Boltzmann factor with the microcanonical distribution, 
but this introduces a fundamental difficulty related to 
the finite (or infinitesimal) support of the latter psf. 

The present work aims at providing a solution to the 
above difficulty by considering an initially microcanoni- 
cal system that evolves in time under a non-Hamiltonian, 
isoenergetic dynamics [2^. As we shall see by explicit 
construction of such a dynamics, the problem associated 
with the microcanonical distribution is overcome, and 
one is able to derive a nonequilibrium identity for entropy 
differences at constant energy, much like Eq. |^ for free 
energy differences at constant temperature [see, in partic- 
ular, Eq. (|19|) ]. In addition, its quasistatic limit reduces 
to an identity that generalizes a well-known thermody- 
namic integration formula for free energies [Eq. (|12(l ]. 
Though admittedly difficult to realize in real experi- 
ments, the isoenergetic processes that underly these re- 
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suits are easily implemented numerically, as discussed 
shortly. 

The remainder of this paper is organized as follows: 
In Sec. ^ the central nonequilibrium identity is intro- 
duced and used to derive some interesting limiting cases 
in Sec. IIIII The importance of these results in the nu- 
merical estimation of entropies and densities of states is 
discussed in Sec. II VI followed by a summary and outlook 
in Sec. |V| 



II. DERIVATION 

To begin the derivation, assume a purely classical 
framework is appropriate, and let T = (x, p) be a point 
in the system phase space, where x, p represent the 
canonically-conjugate position and momenta coordinates 
of all the particles, respectively. The system of interest 
has a Hamiltonian H = H\{pc,p) that depends para- 
metrically on a predefined time-dependent function A(t). 
This function is an external parameter that "switches" 
from A(0) = ^ to A(r) = i? in a given switching time 
T. In order to represent an isoenergetic process, one 
needs to model a suitable "thermostat" that exchanges 
heat with the system so as to precisely counterbalance 
the work done by the external parameter A. One possi- 
ble way of achieving this is to modify Hamilton's equa- 
tions of motion by introducing an artificial "force" field 
F = {F^,Fp) on r as 



(2) 

P = -— + *pli'j: (3) 
and demand the constancy of the system energy through 



dp 



dH 
~dt 



dH 

9A 



A = 0, 



(4) 



where V = d/dT. Using the modified equations of 
motion, one immediately sees that the condition ex- 
pressed by Eq. Q is fulfilled by any vector F satisfying 
■ F — —XdH/dX, i.e. the desired force field is given 
in general by 



-A 



X dH 



(5) 



where X = X(r) is an arbitrary vector field not com- 
pletely perpendicular to \7H. Particular examples of 
X satisfying this condition are X = {0,dH/dp), and 
X = VH itself. The former resembles the "Gaussian" 
thermostats widely used in the literature to model and 
simulate nonequilibrium processes (see e.g. Refs. [Tll[T2l | 
for reviews), and will soon be discussed in more detail. In 
general, however, Eq. |SJl provides a more flexible recipe 
particularly suited for parameter-dependent Hamiltoni- 
ans. 



The time-dependent ensemble density for the above 
non-Hamiltonian dynamics will now be derived. By di- 
rect integration of Liouville's equation in Lagrangian 
form (cf. Eq. (3.3.8) in Ref. |0), one obtains the general 
solution 



where 



p*(ro = po(ro)e-*^*(r' 



AtiTt) = - / dsA(r, 



(6) 



(7) 



is the time average of the "phase space compression fac- 
tor" A(r) = V • r along the trajectory that connects Fq 
to Ft. The notation At (Ft) implies that we are look- 
ing at the above functional of the phase space trajectory 
{F^,}, s : [0, i] as a function of the endpoint Ft, this be- 
ing possible due to the deterministic character of these 
trajectories. For a system evolving under the isoener- 
getic equations of motion Eqs. (01- ©j with F given by 
Eq. lO, one has 



At (Ft) 



ds AV ■ 



X dH 
X-VH dX 



(8) 



so the desired density at r is completely determined by 
Eqs. ©, ^ and the initial condition 



po(ro) 



d[E-HA{T„)] 

nA{E) 



(9) 



which is just the microcanonical distribution. Here S is 
the Dirac deha function, and QxiE) = J dTS[E-HxiT)] 
is the density of states at the external parameter A. 

Consider now the following average over all realizations 
of a switching process that takes X{t) from A to B in t 
units of time: 



= J dTr p,(F,)e^^-(r-) 
S[E - Ha{To)] 



= / dTr 



rtAiE) 



(10) 



In the second line above, Eqs. © and ^ have been 
used. Since the system Hamiltonian is constant along 
any trajectory generated by Eqs. HJ-© [cf. Eq. Q], in 
particular iJ^(Fo) = HB{Tr), it follows from Eq. l(Tn| 
that 



(11) 



where A,- is given by Eq. and AS{E) = 

\n['[lB{E)/^A{E)] is the entropy difference between the 
thermodynamic states A and B. The identity above, 
along with its quasistatic version [Eq. H12I) ]. are the cen- 
tral results of the present paper [see also Eq. H19|) ]. Some 
computational aspects of these results will be discussed 
shortly. 
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III. LIMITING CASES 

Let us now study the quasistatic limit of Eq. i.e. 
the hmiting case where X{t) changes infinitely slowly from 
A to B. In this case, by a slig ht extension of the "adia- 
batic ergodic hypothesis" [Ijl, one expects that the dy- 
namical state of the system (Fs) spends a sufficiently 
long time sweeping a nearly constant-energy, constant-A 
surface E = H\{T s) so that one can invoke the approxi- 
mation 



1 



ds (A) 



_B,A(s) ' 



(quasistatic) 



where the subscripted angle brackets denote a micro- 
canonical ensemble average at constant E and A(s). Ac- 
cordingly, the time average A,- is independent of the spe- 
cific realization of the switching process, and it follows af- 
ter a simple change of integration variables that Eq. 
reduces to the following new thermodynamic integration 
formula: 



ASiE) 



dX(\7 ■ 



X dH 
X~VH~dX 



(12) 



E.X 



Recall that the Jarzynski equality also reduces to a ther- 
modynamic integration formula in the limit of quasistatic 
processes, viz. AF{T) — J^dX {dH/dX)^ 01, a result 
that can be easily derived by standard statistical me- 
chanics. The existence of an arbitrary (see above) vector 
field in Eq. (|12|l is analogous to a generalization of the 
thermodynamic integration formula for AF also derived 
by Jarzynski [T3| . 

An independent proof of Eq. p2|l based on ergodic ma- 
nipulations is now provided as a consistency check. 
Starting from AS{E) = J^dX{dS/dX) and Sx{E) = 
In fl\{E), one has 



''drs[E - Hx{r) 



ndE 
ndE 



dH 
'dX 



da dH 

\VH\ dX ' 



(13) 



where S is the surface E ~ Hx (F) , and da is an infinites- 
imal area element of this surface. Now let X — X(F) 
be a vector field not completely perpendicular to da = 
daVH/\VH\. Then 



da 



— da. ■ 



X 



X-VH' 



so that, by the divergence theorem and a subsequent en- 
ergy differentiation, Eq. (|13|l becomes 



- V. 



dH 



X-VH dX 



which coincides with Eq. (|12() upon integration, q.e.d. 

Although Eqs. (HTJ and lO are exact and can in 
principle be directly utilized in numerical simulations 
(Sec. IIV|I . it is interesting to investigate their limit in 
the case of large systems. Consider first the integrand of 
Eq. JHJ. One has, exactly. 



X dH 
X-VH~dX 



1 dH X 
T~dX ^ ~ 



SJjdH/dX) 
X-S/H ' 



(14) 



where 



r x-vH 

is an inverse "instantaneous temperature" that coincides 
with the inverse of a recently derived expression for 
the microcanonical temperature IT^ upon ensemble- 
averaging. Assume now that, as expected, T is inten- 
sive. Two possibilities of interest arise: (a) A couples to 
a selected number of particles and X thermostats all the 
particles, or (b) A couples to all the particles and X ther- 
mostats any number of particles [2^. Under these condi- 
tions and as the total number of particles is increased, 
the second term in the r.h.s. of Eq. (|14|1 either (a) van- 
ishes or (b) becomes negligible in comparison to the first 
one. For large enough systems, therefore, Eq. © can be 
approximated as 



A. 



T Jo 



at X — , (large systems) 



and the nonequilibrium identity Eq. (|lll) can be rewritten 
as 



(large systems) 



where XdH/dX has been identified with the rate of 
work input due to the external parameter A [see also 
Eq. igj], i.e. 



d M 

dW = — dX. 
dX 



(15) 



By invoking Jensen's inequality (e^) > e^^^ 
above result can also be recast in a form that resem- 
bles the second law of thermodynamics for isoenergetic 
processes (i.e. the Clausius inequality with dQ = —dW), 
namely 



AsiE)>- 




(large systems) (16) 



where, as in the case of the JE, the angle brackets can 
presumably be dropped in the thermodynamic limit 2^. 
It is not difficult to verify, however, that in the case of 
quasistatic processes for sufficiently large systems, one 
obtains from either Eq. H12|l or Eq. H16|l the thermody- 
namic statement 



(large systems, quasistatic) 



4 



where 1/T = (1/T)^ ^ is the inverse equihbrium temper- 
ature of the system, and hence of the reservoir along the 
isoenergetic path. 



IV. UTILITY AND PRACTICAL ASPECTS 

It is natural to inquire about the utility of the above 
results, since they were derived after the admittedly ar- 
tificial equations of motion described by Eqs. 10)- Q- 
Though in physical experiments it is very unlikely that 
one will ever encounter such an isoenergetic set ting, this 
is not the case in numerical simulations (see e.g. [lllll^ ). 
These results should therefore be of greatest interest for 
the computational community, but even in this case one 
might ask why the computation of S{E) or ^1{E) is at all 
relevant, since quantities such as free energies at constant 
temperature/pressure have a more direct connection with 
experiments. 

Although an immediate answer to the above question 
can be found in the context of phase transitions of finite 
or "small" systems I '^i^l briefly describe how 

one can recover important thermodynamic quantities in 
the isothermal ensemble given the knowledge of S{E) or 
D,{E) — e'^'-^-' over a range of energy values. How to 
efficiently obtain entropy and density of states at more 
than a single energy value will be discussed shortly. 

First, notice that the canonical ensemble average at 
temperature T of any energy-dependent observable f{E) 
can be written as 



{fiE))T = 



J^dE'n{E')e-^'/^fiE') 
J^dE'n{E')e'E'/^ 



(17) 



Therefore, knowledge of the density of states over the 
relevant range of energies where Q(E') e~^ 1'^ is concen- 
trated allows one to calculate observables such as average 
potential energy, heat capacities, and free energy differ- 
ences by simple quadrature. In fact, ^(Ei) plays a central 
role in the so-called "flat histogram" Monte Carlo meth- 
ods 1^ , and is used precisely as described above to 
estimate such observables. It is worth mentioning that 
the present nonequilibrium method shares an additional 
feature with flat histogram methods, which have been 
introduced with the goal of overcoming energy barriers 
in finite-temperature simulations. An estimate based on 
Eq. (|ll|l also has the potential of sampling the phase 
space without becoming trapped in energy basins defined 
by barriers much greater than T (recall that fc^ = 1 in 
the present units). This follows from the ability to con- 
struct a switching process that changes the strength of 
the interaction among the particles, as described in the 
next paragraph, so that during the course of the simula- 
tion the dynamics enjoys much lower energy barriers and 
hence greater mobility. 

The results derived in this work can be made highly 
efficient by borrowing a scaling idea originally developed 
for the adiabatic switching method 23]. In the present 



language, the scaling idea allows one to obtain entropy 
differences or ratios of densities of states between the sys- 
tem of interest and a reference system over a wide range 
of energies from the data obtained at a single energy 
value E. In fact, assuming that the reference system is 
an ideal gas, i.e. the parameter-dependent Hamiltonian 
is Hx{x,p) = p^/2 -I- AJ7(x) with A = X{t) switching 
from to 1 in r units of time, by rearranging the expres- 
sion ^1\{E) = Jdxdp S[E — p^/2 — AC/ (x)] one can easily 
derive the following scaling relation: 

SiiE/X) - SoiE) = ASxiE) - In A. (18) 

Here Si{E/X) is the entropy of the system of interest at 
the energy value E/X, So{E) is the (known) ideal gas 
entropy, and the difference ASx{E) = Sx{E) - Sq{E) 
is obtained either directly from Eq. (|12|l with the upper 
integration limit replaced by A, or by recording the inter- 
mediate values of the average in Eq. (|ll|l at the instant t 
corresponding to the desired A — X{t). Note finally that 
the expressions derived herein for AS{E) can equally 
well be used for estimating ratios of density of states, 
ils (£■) /f^A (£-) , without resorting to the traditional his- 
tograms. 

It follows from the above scaling identity that, pro- 
vided the entropy /density of states of the reference sys- 
tem is known with great accuracy, the numeric value of 
S{E) or ri(i?) can be obtained for a wide range of energies 
from a single isoenergetic simulation and, in particular, 
one is able to use Eq. (fTTjl to estimate various thermo- 
dynamic quantities of interest over a range of tempera- 
ture values. In practice, of course, one has to choose the 
range of E judiciously so that the estimates of n{E) are 
reasonably accurate in the range where the integrand of 
Eq. l(T7|) is greatest. 

Several other properties enjoyed by the present results 
will now be discussed. In contrast to the "adiabatic 
switching" method for computing entropy differences , 
the exact result in Eq. Hll() does not suffer from the sys- 
tematic errors associated with a finite switching time |l| ■ 
Moreover, since neither Eq. pi(l nor H12|l relies on a single 
dynamical trajectory, these equations are trivially par- 
allelizable. Isoenergetic equations of motion similar to 
those adopted in this work have been extensively used 
in the literature as means of simulating nonequilibrium 
processes !ll| '24^, and the general recipe provided in 
Eqs. (0)-® and ((Sj should not pose any additional tech- 
nical difficulty. Nonetheless, a considerable operational 
simplification is achieved when X — (0, dH/dp) — (0, p), 
where it is assumed that Hx{x., p) = p"^ /2+Ux{x). In this 
case, the thermostatting mechanism reduces to a sim- 
ple velocity-dependent force Fp = —X{dU/dX)p/p^, and 
Eq. lfTT|) yields exactly 



-J^dW/T 



(19) 



where dW = dtX{dU/dX) is the infinitesimal amount 
of work performed on the system, and 1/T = {Nd — 
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2)/p^ is the inverse instantaneous temperature defined 
in Sec. 11111 {d is the number of spatial dimensions and N 
is the number of particles). The above equation has the 
advantage of involving less abstract and more physically 
sound quantities, while remaining as rigorous as the more 
general case of Eq. (fTT|l . 

V. CONCLUSIONS 

In summary, the present contribution has introduced 
two fundamental nonequilibrium and equilibrium equal- 
ities for entropy differences (or ratios of densities of 
states), Eqs. and (|12|l respectively [see also 

Eq. (|19|l ]. The former is a nontrivial extension of the 
Jarzynski equality [Eq. (Q] to the microcanonical ensem- 
ble, which was made possible through the construction 
of a suitable set of non-Hamiltonian equations of motion. 
This particular extension circumvents the mathematical 
difficulties that one encounters by adapting straightfor- 
ward strategies that have been successful in other en- 
sembles (9i, lOj. The latter is the entropic analog of the 
well-known thermodynamic integration formula for free 
energies. 

The utility of these results was discussed in Sec. IIVI 
where it was remarked that the artificial character of the 



isoenergetic equations of motion described by Eqs. (0)-® 
causes such results to be of greatest interest in numeri- 
cal simulations. In this context, the identities derived in 
this work allow one to compute S{E) (entropy at a given 
energy) or n{E) (density of states) for a wide range of en- 
ergies from a single isoenergetic simulation, as discussed 
in connection with Eq. H18|) . Knowledge of these quan- 
tities per se is of great interest in the study of phase 
transitions of finite or "small" systems , but can 

also be used to recover important observables in other 
ensembles (e.g. isothermal), as evidenced by Eq. (|17|) . 

Since Eqs. pil) and 112|) share in common an arbitrari- 
ness through the vector field X, a question that deserves 
further investigation is whether one can find an optimal 
form for X that maximizes computational efficiency. 
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